#!/usr/bin/perl
#Programmer: Rory Carmichael
#Purpose: Get mean and standard deviation for %Min-Max results for each position of the provided random reverse translations.  Also create a tiny summary file consisting of the mean and standard deviation over all sequences at all positions.
#Usage: rrt-stats.pl <sorted reverse translation file>

#my $rrtf = $ARGV[0];
#open(IF, $rrtf);
# I guess we need to build arrays for each position of each unique sequence.  This is going to be godawful
my $curid;
my $curlength;
my @averages;
my @values;
my @stdevs;
my @allvals;
while (my $line = <>) {
	chomp $line;
	if ($line =~ /^>/) {
		if ($line ne $curid) {
#			print STDERR "doing $curid stats\n";
			if (defined $curid) {
				for(my $i=0;$i<$curlength;$i++) {
					for(my $j=0;$j<$#values+1;$j++) {
						$averages[$i] += $values[$j][$i];
					}
					$averages[$i] /= $#values+1;
				}
				for(my $i=0;$i<$curlength;$i++) {
					for(my $j=0;$j<$#values+1;$j++) {
						$stdevs[$i] += ($values[$j][$i] - $averages[$i]) ** 2;
					}
					$stdevs[$i] /= $#values;
					$stdevs[$i] = sqrt($stdevs[$i]);
				}
				print "$curid avgs\n";
				print join(',', @averages) . "\n";
				print "$curid stdevs\n";
				print join(',', @stdevs) . "\n";
			}
#			print STDERR "initializing array for $curid\n";
			$curid = $line;
			$curlength = 0;
			@values = ();
			@averages = ();
			@stdevs = ();
			my $nextline = <>;
			chomp $nextline;
			my @arr = split(/,/,$nextline);
			$curlength = scalar(@arr);
			push(@values,  [ @arr ]);
			push(@allvals, @arr);
		} else {
#			print STDERR "adding to $curid stats\n";
			my $nextline = <>;
			my @arr = split(/,/,$nextline);
			push(@values, [ @arr ]);
			push(@allvals, @arr);
		}
	} else {
		print STDERR "NOT A TWO LINE SEQUENCE FILE\n";
	}
}
#print STDERR "doing $curid stats\n";
for(my $i=0;$i<$curlength;$i++) {
	for(my $j=0;$j<$#values+1;$j++) {
		$averages[$i] += $values[$j][$i];
	}
	$averages[$i] /= $#values+1;
}
for(my $i=0;$i<$curlength;$i++) {
	for(my $j=0;$j<$#values+1;$j++) {
		$stdevs[$i] += ($values[$j][$i] - $averages[$i]) ** 2;
	}
	$stdevs[$i] /= $#values;
	$stdevs[$i] = sqrt($stdevs[$i]);
}
#print "$curlength was curlength\n";
#print $#values+1 . " sequences evaluated\n";
print "$curid avgs\n";
print join(',', @averages) . "\n";
print "$curid stdevs\n";
print join(',', @stdevs) . "\n";
my $allavg = 0;
my $allstdev = 0;
foreach my $val (@allvals) {
	$allavg += $val;
}
$allavg = $allavg/($#allvals+1);
print STDERR "Total Average: $allavg\n";
foreach my $val (@allvals) {
	$allstdev += ($val - $allavg) ** 2;
}
$allstdev /= $#allvals;
$allstdev = sqrt($allstdev);
print STDERR "Total STDEV: $allstdev\n";
#print STDERR "Suggested Cutoff: " . $allavg - 7*($allstdev) . "\n";
